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ABSTRACT 

We study the process of cosmic reionization and estimate the ionizing background 
in the intergalactic medium (IGM) using the Lyman series absorption in the spectra 
of the four quasars at 5.7 < z < 6.3 discovered by the Sloan Digital Sky Survey. 
We derive the redshift evolution of the ionizing background at high redshifts, using 
both semi-analytic techniques and cosmological simulations to model the density 
fluctuations in the IGM. The existence of the complete Lya Gunn-Peterson trough in 
the spectrum of the z = 6.28 quasar SDSS 1030+0524 indicates a photoionization rate 
(r_i2 in units of 10~^'^s~^) at z ~ 6 lower than 0.08, at least a factor of 6 smaller than 
the value at z ~ 3. The Ly/3 and Ly7 Gunn-Peterson troughs give an even stronger 
limit r_i2 <^ 0.02 due to their smaller oscillator strengths, indicating that the ionizing 
background in the IGM at z ~ 6 is more than 20 times lower than that at z ~ 3. 
Meanwhile, the volume-averaged neutral hydrogen fraction increases from 10~^ at 
2: ~ 3 to > 10~^ at z ~ 6. At this redshift, the mass-averaged neutral hydrogen fraction 
is larger than 1%; the mildly overdense regions {5 > 3) are still mostly neutral and the 
comoving mean free path of ionizing photons is shorter than 8 Mpc. Comparison with 
simulations of cosmological reionization shows that the observed properties of the IGM 
at z ~ 6 are typical of those in the era at the end of the overlap stage of reionization 
when the individual HII regions merge. Thus, z ~ 6 marks the end of the reionization 
epoch. The redshift of reionization constrains the small scale power of the mass density 
fluctuations and the star forming efficiency of the first generation of objects. 
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1. Introduction 

After the recombination epoch at z ~ 1500, the universe remained mostly neutral, until 
the first generation of stars and quasars ionized the intergalactic medium (IGM) and ended the 



cosmic "dark ages" (e.g. Rees 199q ). Popular cosmological models predict this reionization to 



have occurred at redshift between 6 and 20 (e.g., Gnedin &: Ostriker 19971 , Chiu & Ostriker 



2000, Gnedin 2000, Ciradi et al. 2001, Razoumov et al. 2001 and references therein). When 



and how the universe reionized is one of the fundamental questions of modern cosmology (for a 



complete review, see Barkana & Loeb 2001 and Loeb & Barkana 2001). Observations of absorption 



spectra of luminous sources at high redshifts approaching the reionization epoch provide the best 
observational probes of reionization to date. 



Recent discoveries of luminous quasars at z > 5.7 (|Fan et al. 2000| , [Fan et al. 200Tc|) usmg 



imaging data from the Sloan Digital Sky Survey (SDSS, |York et al. 2000 ) enable us to study 



the state of the intergalactic medium (IGM) at redshift up to 6. Observations of Fan et al. 
(2001c, hereafter Paper I), Becker et al. (2001, hereafter Paper II) and Pentericci et al. (2001, 
hereafter Paper III) show that the Lya absorption due to neutral hydrogen in the IGM increases 
dramatically toward high redshifts. In particular, Keck and VLT spectroscopy of high redshift 
quasars (in Papers II and III, respectively) show the first observation of a complete Gunn-Peterson 
(GP) trough (Shklovsky 1964, Scheuer 1965, Gunn & Peterson 1965) in the spectrum of the 
z = 6.28 quasar, SDSSp J103027.10+052455.0 (SDSS 1030+0524 for brevity), where no flux is 
detected over a 300A region immediately blueward of the Lya emission line. The flux decrement 
between the red and blue sides of the Lya emission line is larger than 150, corresponding to an 
effective optical depth to Lya photons (reg) larger than 5 at z^^^ ~ 6.05. The existence of a Ly/3 
GP trough in the Keck spectrum imposes an even stronger limit on the effective equivalent Lya 
optical depth (rgfr > 20, Paper II). pjorgovski et al. (2001) also observed a particularly dark 



region of length ~ 5 Mpc at z ~ 5.4, along the line of sight to the z ~ 5.8 quasar SDSS 1044-0125. 

This is a theoretical companion paper to Papers I - III. In this paper, we use the absorption 
measurements presented in the previous papers to calculate the evolution of the ionizing 
background and the neutral hydrogen fraction in the IGM, and constrain the epoch of reionization. 
The outline of this paper is as follows. In §2, we summarize the observed evolution of the average 
Lya absorption along the lines of the sight to the sample of quasars in Paper I based on the spectra 
presented in Papers II and III, and discuss the errors on these measurements. In §3, we calculate 
the ionizing background at different redshifts using the average Lya absorption, taking into 
account the inhomogeneity of the IGM using a model for the overdensity distribution. We derive 
stronger constraints on the ionizing background in the GP trough region of SDSS 1030+0524 at 
-2abs ~ 6 using the Ly/3 and Ly7 absorption troughs. In §4, we use N-body simulations to model 
the high-redshift Lya forest spectra, and compare the statistical properties of the simulated and 
observed spectra. In §5, we calculate the redshift evolution of the volume and mass weighted 
neutral hydrogen fractions in the IGM, and show how the reionization penetrates into progressively 
overdense regions at lower redshifts as the mean free path of ionizing photons increases. These 
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calculations are used to constrain the epoch of reionization. We discuss the prediction of structure 
formation models for the reionization epoch, and outline future observations needed to probe deep 
into and beyond the reionization epoch in §6. 

Throughout the paper, we adopt a LCDM cosmology with = 0.35, A = 0.65, h = 0.65, 
and Qi)h^ = 0.02 unless otherwise noted. 



2. Evolution of Neutral Hydrogen Absorption 

The Gunn- Peterson (1965) optical depth to Lya photons is 

TGP = faKH'^{z)nu_i, (1) 

nieC 

where fa is the oscillator strength of the Lya transition, = 1216A, H{z) is the Hubble 
constant at redshift z, and nni is the density of neutral hydrogen in the IGM. At high redshifts, 
H{z) oc hQm'^{l + z)^^'^, and the GP optical depth for a uniformly distributed IGM can be 
re- written as: 

In reality, the IGM is highly inhomogeneous, and regions with different densities will be in different 
ionization states and have different GP optical depths (§3). The IGM density distribution needs 
to be taken into account when estimating the averaged neutral fraction from the observations (§5). 

Figure 1 shows the observed evolution of the average Lya absorption in the high-redshift 
quasar spectra as a function of redshift, both in terms of the transmitted flux ratio T, defined as 
the ratio of observed and unabsorbed continuum fluxes (see Paper I), 

T(Zabs) = (/f ) , (3) 

and the effective GP optical depth, defined as 

Taf. = -ln(r). (4) 



The measurements shown in Figure 1 at z^hs < 4.5 are taken from McDonald &: Miralda- 



Escude (2001) . The measurements at 4.8 < Zabs < 5.6 are based on the spectra of Paper II and 
III. In this Figure, we have averaged the measurements along four lines of sights into four different 
redshift bins, Zabs = [4.8 - 5.2], [5.2 - 5.6], [5.6 - 5.95] and [5.95 - 6.15], respectively, and include 
the systematic errors from continuum shape and sample variance (see below). We compute these 
quantities using only the data between Zabs < -^em — 0.1 and 1 + Zabs > 1040 x (1 + 2;cm)/1216, to 
avoid contamination from the quasar Lya emission line and the proximity effect on the red side, 
and from the quasar Ly/3 +OVI emission lines on the blue side. The value in the highest redshift 
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bin, z = [5.95 - 6.15], T = 0.004 ± 0.003, is consistent with zero transmitted flux. Figure 1(b) 
shows the l-cr lower Hmit on the optical depth r^p > 5.1 from the Lya measurement. 

The error on the transmitted flux ratio includes contributions from three sources: (1) photon 

noise in the spectrum, (T(photon), which is measured directly from the data. The typical (7(photon) 
is of the order 2 x 10^^. (2) systematic uncertainties due to the extrapolation of the intrinsic 
quasar continuum, c7(con). In Paper I-III, we assumed a power law continuum fi, oc where 
a = 0.5, normalized at rest frame wavelength 1280A. Therefore, 



log 



1216(1 +Zabs) 



1280(1 + Zc 



(j(a) X r. (5) 



For a{a) ~ 0.4 (Fan et al. 2001a), the typical relative error is acon/T ^ 0.05. (3) sample variance, 
(j(sample). This quantity can be roughly estimated using a model in which the Lya clouds are 
distributed randomly in the IGM (Zuo 1992). We calculate (T(sample) using Equations (8) and (9) 
of Zuo (1992), assuming that the column density distribution of the Lya clouds follows a power law 
N{Nhj) oc , with /3 = —1.5, and the Doppler width of the Lya forest line is 6 = 30 km s~^. In 
a rcdshift window Sz = 0.4, cr(samplc) 0.025 for T ~ 0.10 while for Sz = 0.2, (7(sample) ~ 0.002 
for T ~ 0.006. Therefore, for all the measurements except at the highest redshift (^abs ~ 6), the 
error term is dominated by sample variance. This conclusion is consistent with the scatter of T 
measurements from different lines of sight at the same redshift, which is larger than that expected 
from photon noise alone. For T <C 1, the contribution from photon noise becomes important. The 
error bars in Figure 1 represent the sum of the three error terms, added in quadrature. 



3. Evolution of the Ionizing Background 

The Lya forest arises from absorption of Lya photons by neutral hydrogen gas in the 
intergalactic medium (eq. 2). Assuming local photoionization-recombination equilibrium, we have: 

nniF = nminea{T), (6) 

where nni, ?^hii and arc the local densities of neutral, ionized hydrogen and electrons in the 
IGM, respectively, T is the photoionization rate, and a(T) is the recombination coefficient at 
temperature T (Abel et al. 1997), 

a(T) = 4.2 X 10-^3(T/10^K)-°-'^cm3s-^ (7) 

The photoionization rate T is related to the ionizing background flux J^, by, 

r = 4./^.„<i., (8) 

where the integral includes the ionizing photons from the HI Lyman Limit to the Hell Lyman 
Limit assuming that the He in the IGM is singly ionized, and ai, is the HI cross section of 
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ionizing photons, a,y ^ u ^. For a power law spectrum of the ionizing background dominated 
by quasar, Ji, oc z^~^-^ ( Madau, Haardt, &: Rees 1999| ), we find r_i2 = 2.5J_2i, where r_i2 



is the photoionization rate in units of 10~^^ s~^, and J_2i is the ionizing flux at the Lyman 
Limit in units of 10~^^ erg cm~^ s~^ Hz~^ sr~^. If the ionizing background at high-redshifts is 
dominated by massive stars (e.g. Fan et al. 2001c) with oc i'' ( Parkana &: Loeb 200 1|) , we find 
r_i2 = 1.5J-_2i- 

If the IGM is mostly ionized by a uniform ionizing background, the evolution of the optical 
depth can be expressed as (Weinberg et al. 1997), 

°^ fWiz) hnW^ • 

The Lya absorption increases rapidly with increasing redshift even if the ionizing background 



remains constant with redshift. McDonald &: Miralda-Escude (2001) estimate the ionizing 



background at z < 5.2 by comparing the observed transmitted flux ratio to that of artificial Lya 
forest spectra created from cosmological simulations. In this section, we develop a semi-analytic 
model which uses an empirical model for the density distribution derived from hydrodynamical 
simulations similar to that of McDonald &i Miralda-Escude (2001) , and use it to estimate the 



ionizing background from the Lya observations (§3.1). Then, we use the Ly/3 and Ly7 forest GP 
trough measurements in the spectrum of SDSS 1030-1-0524 to put stronger constraints on the 
background at z ~ 6 (§3.2). 



3.1. Semi-analytic Modelling of Lya Absorption 



The Lya forest arises from low-density gas in the IGM that is in approximate thermal 
equilibrium between photoionization heating by the UV background and adiabatic cooling due to 
Hubble expansion (see e.g., Bi 1993 ; pen et al. 1994| ; [Zhang, Anninos, &: Norman 1995 ; Hernquist 
let al. 19"96t Hui fc Gnedin 19971) . The neutral hydrogen fraction and therefore the GP optical 
depth depends on the local density of the IGM. We define the fractional density of the IGM as 
A = p/{p) (=6 + 1, where 5, the density contrast, is the departure of the local density from the 
mean density, in units of the mean density). For a region of IGM with density A 



"^^^ " hnA,z)noj — ^ • 



(10) 



The dependence on A^ arises because r oc nm, which is proportional to n^jj (eq.6), and 
proportional to A^ for a highly ionized IGM. The temperature of the IGM is determined 
by photoionization-recombination equilibrium, which leads to a power-law relation between 
temperature and density, of the form T = TqA'^, with Tq ~ 1 — 2 x 10^ K, and 7 ~ — 1 (e.g. Hui 



fe Gnedin 1997). Following McDonald fc Miralda-Escude (2001) , we assume an uniform ionizing 
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background, and 7 = 0. Thus, 



/ ^ + zV-^ ( 0.05 , 
r(A) = TO -— A^. (11) 



7 ; Vr_i2(z) 

We determine tq below. The mean transmitted flux ratio T can be calculated as 

POO 

T={e-^)= / e-^(%(A)(i(A), (12) 

JO 

where p(A) is the distribution function of the density of the IGM. For an inhomogeneous IGM, 
using the definition of effective optical depth in Eq.(4) and (12), we have T = e~'^°^ = (e~^) > e~^'^\ 
thus Teff < (r) . 

We calculate T using a parametric form for the volume-weighted density distribution function 
p{A) (Miralda-Escude, Haehnelt, & Rees, 2000): 



p(A) = j4exp 



(A-2/3 _ Cq)2 

2(2(5o/3)2 



A-^ (13) 



where 6q = 7.61/(1 + z), and (3, Cq and A are numerical constants given in Table 1 of Miralda- 
Escude et al. (2000) at several redshifts. This distribution is derived assuming that the initial 
density fluctuations form a Gaussian random field, the gas in voids is expanding at a constant 
velocity, and the density field is smoothed on the Jeans scale of the photoionized gas. It reproduces 
well the density distribution of the photoionized gas in the LCDM hydrodynamic simulations of 
Miralda-Escude et al. (1996). In the following calculations, we linearly interpolate f3 in redshift, 
and set the constants Cq and A so that both the total volume and mass are normalized to unity. 
Finally, we find the numerical constant in Eq.(ll) to be tq = 82 by requiring the resulting T 



to match the measurement of McDonald fc Miralda-Escude (2001) at z = 4.5. Note that the 



quantity we are really constraining with Lya absorption is the combination TQm h~ Tq 
(equations 10 and 11). There are still considerable uncertainties in the values of ^Im, h and 
To, all of which translate to uncertainties in the normalization of ri2. However, in the era after 
reionization, the temperature of the IGM (Tq) is not a strong function of redshift. Therefore, the 
redshift evolution rate of derived here is robust to uncertainties in the normalization. Note 
that [McDonald fc Miralda-Escude (200ll assume that Tq = 2 x lO^K, n^h'^ = 0.02, h = 0.65 and 



= 0.4. The Tuiz) would be ~ 40% higher for Tq = lO'^K (e.g., [Schaye et al. 2000| ). [McDonald 



fc Miralda-Escude (2001) show that the estimate of is quite insensitive to the slope of the 



A — T relation (i.e., the value of 7). 

Also note that the effective optical depth here is calculated by integrating in real space rather 
than in redshift space. We compared the results from N-body simulations using redshift and real 
space coordinates (§4) and found that the effect of ignoring peculiar velocities on the resulting 
optical depth is negligible. 



Figure 2 shows the evolution of the photoionization rate r_i2 based on the observations 
presented in Figure 1, including the measurements at z < 4.5 from McDonald & Miralda-Escude 
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(2001). The ionizing background declines toward higher redshift, from r_i2 ^ 0.5 at Zabs ~ 3, to 
r_i2 ~ 0.2 at Zabs ~ 5, and to r_i2 ~ 0.12 at Zabs ~ 5.8. Note that F remains roughly constant 



between z ~ 4.4 — 5.5, albeit with rather large error bars on individual measurements (see also Cen 
fc McDonald 2001[ ). In the GP trough region in the highest redshift object, SDSS 1030+0524, the 
upper limit on the transmitted flux ratio, T < 0.006, yields a much more stringent upper limit on 
the ionizing background at ^abs ~ 6.05: r_i2 < 0.08, a factor of 6 smaller than the measured value 
at Zahs ~ 3. This dramatic decline in the ionizing background indicates a much larger neutral 
fraction of the IGM and a much shorter mean free path of the ionizing photons at Zabs = 6 that at 
^abs = 3 (§5). 



3.2. Constraints from Ly/? and Lyj Troughs in SDSS 1030+0524 

As pointed out in Paper II, because of its much smaller oscillator strength, and the fact that 
the flux scales exponentially with oscillator strength, one can put a stronger limit on the neutral 
hydrogen fraction in the GP trough from the Ly/3 transition than from the Lya transition, even 
after allowing for the foreground Lya absorption at lower redshifts. In this subsection, we put a 
much stronger limit on r_i2 based on the Ly/J measurement of Paper II. We also attempt to use 
the Ly7 trough to set a limit on r_i2, although it suffers from larger uncertainties in the modelling 
due to the presence of both foreground Ly(3 and Lya absorption, and it is contaminated by Ly6 
absorption. 

The Ly/3 absorption at redshift Z/3 overlaps with the foreground Lya absorption at the redshift 

Za = iXa/Xp){l + Zp)-l. (14) 

where Xa,/3 are the rest-frame wavelengths of the Lya and Ly/3 transitions, respectively. The 

observed total optical depth at the corresponding wavelength is rtotai = Taiza) + Tp{zp). For a 
fractional density A, 

r^(A) = ifp/UMA) = r,(A)/5.27, (15) 

where fa and fp are the oscillator strengths of Lya and Ly/3 transitions, respectively, and ra(A) 
can be calculated using Eq.(ll). Therefore, the observed transmitted flux ratio in the presence of 
both Lya and Ly/? absorption can be written as, 

Ta+p = (e---->) =/^J^^e-[-"(^-^")+-M^^'^/')p(z„,A,)p(z^,A^)dA,(iA^ 

= [/a. e-<^(--^<^)p(z,, Aa)dAa] [/a^ e-^(^^)pizf,, Ap)dAp] (16) 

= Ta{Za)Tf3{zi3), 

where Ta{za) is the transmitted flux of the foreground Lya absorption at z^ defined in Eq.(13), 
and Tji{zfj) is the transmitted flux due solely to the Ly/3 absorption. 

In Paper II, we found = -0.002 ± 0.020 in the Ly/3 GP trough of SDSS 1030+0524 at 
-Zabs = 6, after correcting for the foreground Lya absorption (7^ = 0.12) at Zabs ~ 5.0. Using 
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Eqs.(15) and (16) above, we find r_i2 < 0.025 from the Ly/? trough, an upper limit a factor of ~ 3 
lower than that estimated from the Lya trough. 

The oscillator strength of the hyy transition is a factor of 14.3 smaller than that of Lya 
and a factor of 2.7 lower than that of Ly/3 ; therefore in principle, Ly7 can be used to put an 
even stronger constraint on the equivalent Ly a optical depth. In practice, the modelling of Lyj 
absorption is complicated due to three reasons: 

1. Foreground Lya absorption at Za ~ 4.6. We assume 7^ ~ 0.20, following the values in Figure 
1. 

2. Foreground Ly/3 absorption at ~ 5.7. From Figure 2, we find T^i2{z = 5.7) ~ 0.12, and 
the Ly/3 transmission 7^ = 0.33 using Equations (13) and (15). 

3. Contamination from Ly5 absorption. 

In the redshift range of the GP trough 5.95 < Zabs < 6.16, Ly7 absorption covers a wavelength 
range 6759A < A < 6963A. The hy6 absorption covers a wavelength range 66OOA < A < 68OOA, 
partially overlapping with the Ly7 trough, and the hy6 emission line is at A = 6914A. From the 
Keck/ESI spectrum presented in Paper II, we find that T = 0.0043 it 0.0033 in the wavelength 
range 6759A < A < 6963A, and T = 0.0028 it 0.0041 in the more restricted wavelength range 
6800A< A < 6963A. The latter wavelength range is redward of the Ly6 GP trough. Although part 
of it is still bluer than the Ly6 emission, the amount of LyJ absorption is strongly suppressed due 
to the proximity effect from the quasar. In the wavelength range of 6914A < A < 6963A, totally 
redward of Ly6 emission, we find T = —0.0022 it 0.0087. 

From a combination of these transmitted flux ratio measurements, we adopt the limit 
Ta+js+'y < 0.007 at 1-0" in the Ly7 GP trough region. Correcting for the foreground Lya and Ly/3 
absorptions, we find the transmission solely due to Ly7 to be < 0.10. This transmitted flux 
ratio corresponds to a photoionization rate of r_i2 < 0.020, comparable to the limit derived from 
the Ly/3 GP trough, and representing a factor of > 20 decline from the ionizing background at 
z ~ 3. The rapid decrease of F at z > 5.7 indicates a possible sharp transition in the ionization 
state of the IGM (§5). 



4. Simulating the Lya forest 

In this section, we use an N-body simulation to simulate the absorption spectra of high- 
redshift quasars in the Lya region, taking into account the resolution and noise properties of 
the observations. We fix the photoionization rate in the simulation to reproduce the observed 
mean transmitted flux at each redshift, and compare the statistical properties of the simulated 
spectra including the probability distribution function (PDF) and the threshold crossing statistics 
(Miralda-Escude et al. 1996), with those of the observed spectra. 
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We carry out the simulations in a LCDM model, with 0^ = 0.3, JIa = 0.7, h = 0.65, 
Qbh^ = 0.02, and fSm = 0.9. Here, crgm is the rms density fluctuation in 8h ^ Mpc spheres, chosen 
here to reproduce the observed abundance of clusters at z = ( [White, Efstathiou, fc Frenk 1993 ; 



Eke, Cole, fc Frenk 1996| ). We evolve the density fluctuations using a particle- mesh (PM) N-body 



code that is described in detail in Hennawi et al. (2001). This code uses a staggered mesh to 
compute forces on particles ( Melott 1986 ; Park 1990|) , and uses the leapfrog scheme described 



in Quinn et al. (1997) to integrate the equations of motion. The periodic simulation cube is 
L = 25/i~^Mpc on a side, and uses Np = 256^ particles and an Nm = 512^ force mesh. We set 
up initial density fluctuations using the power spectrum parameterization of Efstathiou, Bond & 
White (1992) with the shape parameter T = 0.2, and assign initial displacements and velocities to 
the particles using the Zeldovich approximation. We evolve from redshift z = 49 — > 5.45, taking 
55 equal steps in the expansion scale factor. 

We use the TIPSYQ package to create artificial Lya absorption spectra along 400 random 
lines of sight through the simulation cube. TIPSY calculates the local mass density (pm) at the 
location of each dark matter particle using a cubic spline smoothing kernel (Hernquist & Katz 
1989) enclosing 32 neighbors, and assigns it a temperature T = TqA"'^. We assumed Tq = 10^ K 
while creating the spectra, although the values of we quote in this paper have been rescaled to 
To = 2 X 10^ K. Spectra along various lines of sight are then created from this particle distribution 
using the algorithm described by Hernquist et al. (1996). We smooth the simulated spectra with 
a Gaussian filter of smoothing length = 28 km s^^ (corresponding to a Full Width at Half 
Maximum of 66 km s~^), and bin them in pixels of width 35 km s~^, approximately twice the 
resolution of the Keck spectra. We then add noise to the simulated Lya forest spectra. The noise 
in the Keck spectra is dominated by the background sky and is a strong function of wavelength 
around 8500A. In order to correctly simulate this structure in the noise, we replace the fluxes 
(F) in each simulated spectrum by the values F — > F + o"obsG'(l), where Cobs is the noise array 
in contiguous pixels starting from a random location in the Keck spectrum binned to about 
35 km s~^, and G(l) is a Gaussian random deviate with zero mean and unit variance. We scale all 
the optical depths at each redshift so that the mean transmitted flux ratio computed from all the 
400 simulated spectra match the observed values. 

We compute the Lya forest spectra at Zabs = 5.5,5.7 and 6.0, and compare them with the 
observed Keck spectra of SDSS 1044-0130, SDSS 1306+0356 and SDSS 1030+0524 at the same 
wavelength ranges. Figure 3 shows the Lya forest at 5.4 < z^hs < 5.6 and 5.9 < Zabs < 6.1 from the 
Keck spectra, and the artificial Lya absorption spectra along eight random lines of sight through 
the simulation. Since the simulation cube is only 25h~^ Mpc on a side, each simulated spectrum is 
of length 3506 km s~^, 3558 km s~^ and 3635 km s~^ at Zabs = 5.5, 5.7 and 6, respectively. The 
simulated spectra are qualitatively similar to the Lya forest region in the Keck spectra. Although 
they have the same mean transmitted flux by construction, the fluctuations in the simulated 



^see http://www-hpcc.astro.washington.edu/tools/tipsy/tipsy.html 
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spectra arise from fluctuations in the mass density field in a LCDM model (except at Zabs = 6 
where both the observed and the simulated spectra are consistent with noise). The fact that the 
fluctuations in the simulated spectra are qualitatively similar to those in the Keck spectra is an 
indication that our physical description of the Lya forest as arising from photoionized gas in a 
low-density IGM in a LCDM model is reasonable. 

Figure 4 shows the probability distribution function (PDF) of the transmitted flux ratio 
( Jenkins &: Ostriker 1991 ; Miralda-Escude et al. 199"^ ; Ranch et al. 1997 ), computed using the 
400 artificial spectra. The solid points in the three panels show the fiux PDF measured from the 
Keck spectra of SDSS 1044-0125, SDSS 1306+0356 and SDSS 1030+0524 in a velocity range of 
8500 km s~^ around (-Zabs) = 5.5, 5.7 and 6.0, respectively. The solid line shows the average fiux 
PDF from the simulated spectra. Since the simulated spectra are shorter in length than the data, 
we compute the cosmic variance in the flux PDF measured from a single Keck spectrum of length 
8500 km s~^ by combining together 2.5 artiflcial spectra end to end so that the resulting artiflcial 
spectrum has roughly the same length as the Lya forest region in the Keck spectrum. The flux 
PDF of the simulated spectra is consistent with that computed from the Keck spectra, in all 
three redshift bins. In order to match the observed mean transmitted flux ratio, we flnd that the 
ionizing background at z ~ 5.5 and 5.7 is r_i2 = 0.21, a factor of 2.5 lower than the value at 
z ~ 3, and consistent with the result from §3.1. Note that the shape of this distribution is related 
to the shape of the power spectrum, and is a prediction from the LCDM model. For z ~ 6.0, the 
simulated spectrum reproduces the observed CP trough with r_i2 < 0.08, again consistent with 
the limit derived in §3.1 using a semi-analytic model. At this redshift, since there is no detected 
flux in the Keck spectrum, the observed and simulated PDFs merely reflect the distribution of 
noise in the spectrum. 

The three panels in Figure 5 show the number of flux threshold crossings per unit redshift 
(the threshold crossing statistic; Miralda-Escude et al. 199"^ ) in the simulated spectra (solid line) 
and the Keck spectra (points) in three redshift bins. We show this quantity as a function of the 



fraction of pixels in the spectrum that are less than the threshold flux (Miralda-Escude et al. 



1996| ; [Weinberg et al. 1998|) . This statistic is analogous to the "genus curve" used to characterize 



the topology of the 3-dimensional galaxy distribution and has the advantage that it is insensitive 
to the exact relation between the relative distribution of dark matter and baryons (the "biasing" 
of the Lya forest), as long as this relation is monotonic. We compute the cosmic variance in this 
statistic in a similar manner to the flux PDF. The simulated spectra reproduce this statistic of 
the observations very well, at all three redshifts. Therefore this modeling of the Lya forest in a 
LCDM model can reproduce the topology of the dark matter density fluctuations that are probed 
by the observed Lya forest in these quasar spectra. 
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5. Constraining the Epoch of Reionization 

In this section, we expand the model described in §3.1 to calculate from observations the 
evolution of three quantities that reflect the ionization state of the IGM: (1) the neutral fraction of 
the IGM, taking into account the inhomogeneity of the IGM, (2) the critical density - the density 
lower than which the gas is completely ionized, and (3) the mean free path of ionizing photons in 
the IGM. Using these calculations, we show that as reionization progresses toward lower redshifts, 
the neutral fraction decreases, regions of progressively higher densities become ionized, and the 
mean free path of ionizing photons increases. We compare these three quantities with the results 
from a simulation of the cosmological reionization process by Gnedin (2000), to constrain the 
likely epoch of reionization (see also Gnedin 2001| ). 



We can use Eq.(6) to calculate the local neutral hydrogen fraction [/hi(A)] in a region of 
density A, 

f (A) = nm{/^) ^ nHi(A) 

^"^^ > nHi(A) + nHii(A) (nH)A ' ^ 
where (nn) is the average nn over the universe. We define 

A = 1.16{nn)a{T)/T, (18) 

where the coefficient 1.16 comes from the helium contribution to the electron density assuming a 
helium fraction Y = 0.24 (Weinberg et al. 1997). Then, it can be shown that 

^ ' 1 + V1 + 4AA ^ ^ 

Note that for a highly ionized medium, /hi ~ ^A. We calculate the redshift evolution of the 
HI fraction using the density distribution of Eq.(13), and the results on T in Figure 2. Note 
that /hi is directly related to A ~ a(T)/T. Therefore, our estimate of /hi is insensitive to our 
assumed value for the IGM temperature T, since our estimate of T is itself proportional to the 
assumed a{T). Figure 6 shows both the volume-averaged neutral hydrogen fraction /^j (dashed 
line) and the mass-averaged quantity /^ (solid line). The volume-averaged HI fraction increases 
from 1.7 X 10~^ at z ~ 3.9 to 1.3 x 10~^ at z = 5.8 ~ an increase by a factor of ~ 8. In the GP 
trough region at z ~ 6.05, the upper limit on the ionization rate F from the Lya , Ly/3 and Ly7 
transitions yields a lower limit on the neutral fraction, ^ 10^^, and /^j ^ 10~^, an increase by 
almost two orders of magnitude from z ~ 4. The mass-averaged value is larger because most HI is 
concentrated in dense clumps in the IGM which occupy a small volume. 

The neutral hydrogen fraction depends on the detailed density distribution of the IGM. As 
pointed out in §3.1, t^s < (t). Therefore, simply using Eq.(2), which assumes a homogeneous 
IGM, will severely underestimate the neutral fraction. For example, at z = 5.8, Tqp = 3.7 (Figure 
1), and the neutral fraction derived using Eq.(2) is 5.6 x 10~^; a factor of ~ 2.3 and ~ 70 lower 
than the volume and mass-averaged neutral fractions, respectively, from Eq.(19). In a clumpy 
IGM, the high-density regions quickly become opaque to all Lya photons, while almost all the 
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transmitted flux comes through the deepest voids in the IGM. We iUustrate this point in Figure 
7, which shows the cumulative distributions of transmitted flux as a function of density at z=5.3. 
At this redshift, while Tcff ~ 2.5, the region with the minimum absorption in the spectrum of 
SDSS 1044-0125 (Fan et al. 2000) has r < 0.5. It represents the deepest void in the IGM, with 
A ~ 0.15 (from Eq. 11), and has a neutral fraction of ~ 1.2 x 10~^, a factor of seven lower than 
the volume-averaged neutral fraction at this redshift. It is a factor of ~ 50 lower than the value 
one would get by directly applying Eq.(2) and assuming a uniform IGM, as r scales as (Eq. 
11). 

In all the calculations above, we have assumed that the ionizing background (and hence 
r) is uniformly distributed in the IGM. Figure 7 shows that the F derived from the observed 
transmitted flux ratio is characteristic of the background level in low density regions of the IGM, 
since the high density regions contribute little to the observed flux. The observed transmitted flux 
[Eq.(ll)] is essentially a volume-averaged measurement. In a similar sense, our calculation of the 
volume-averaged neutral fraction should be regarded as more reliable than the mass-average, as 
the former quantity is less affected by the the local ionization state in the high density regions. 

Miralda-Escude et al. (2000) present a semi-analytic model describing the evolution of 
reionization in an inhomogeneous universe; reionization starts in voids and then gradually 

penetrates deeper into ovcrdcnse regions. This process is defined by the redshift evolution of the 
critical density parameter Aj. At any redshift, regions with A > Aj remain neutral, while regions 
with A < Aj are mostly ionized. Here we follow Gncdin (2000) by defining Aj as the critical 
density above which 95% of the neutral gas lies. Figure 8 shows the evolution of Aj based on 
the neutral fraction calculated above. The quantity Aj ~ 40 at z ~ 4, and it quickly decreases 
to Aj < 4 in the GP trough region at 2 ~ 6. Using Eq.(8) of Miralda-Escude et al. (2000), we 
calculate the mean free path of the ionizing photons, Aj = Ao[l — Fy{Ai)]~'^/^ , where Fv(Aj) is the 
fraction of volume with A < Aj, and XqH = 60 km s^^ (Miralda-Escude ct al. 2000). Figure 9 
shows the evolution of the comoving mean free path with redshift. It decreases from ~ 80h~^ Mpc 
at ^ ~ 4 to smaller than ~ 8h~^ Mpc in the GP trough region at ^ ~ 6. 

Figures 6, 8 and 9 clearly demonstrate that at z ~ 6, the universe is much more neutral 
than at lower redshift (z ~ 3 — 5). At z ~ 6, the presence of the GP trough shows that more 
than 1% of the protons are in the form of HI, the moderately overdense regions (A ^ 4) are still 
mostly neutral, and the mean free path of ionizing photons is several Mpc in comoving distance 
— a distance comparable to the correlation length of the galaxy distribution at z = 0. The last 
remaining HI in the IGM is being ionized at this epoch. In this sense, the first appearance of the 
complete GP trough at z ~ 6, marks the end of the reionization epoch. 

Cosmological reionization is a complicated process. Cosmological simulations that include 
gas dynamics, star formation processes, atomic and molecular physics and radiative transfer have 
recently begun to to shed light on when and how this process happened in popular cosmological 
models. Gnedin (2000) divides the reionization process into three phases: the pre-overlap stage in 
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which individual HII regions begin to grow; the overlap stage in which these HII regions merge 
and the ionizing background in the IGM rises by a large factor; and the post-overlap stage in 
which the remaining HI in the high density regions is gradually ionized. The results from the 
simulation of Gnedin (2000) are shown in Figures 6, 8 and 9. In the simulation, the overlap stage 
was assumed to be at 2; ~ 7. At this redshift, the universe goes through a near phase-transition: 
the ionizing background and the mean free path of photons increases and the neutral hydrogen 
fraction decreases dramatically over a narrow redshift range. The reionization redshift can be 
defined in a number of ways: e.g., the redshift at which the average neutral hydrogen fraction 
is ~ 50%; the redshift at which the HII filling factor is of order unity; the redshift at which an 
average region of the IGM is ionized (Aj ~ 1), etc. However, as shown in Gnedin (2000), Chiu &; 
Ostriker (2000) and Razoumov et al. (2001), the overlapping stage occurs over a small range in 
redshift, with the neutral fraction changing from 90% to < 1% in Az < 2 (corresponding to a time 
interval of 200 million years at z ~ 7), while r_i2 increases from ~ 10~^ to ~ 10"^. 

Comparing the predictions from the simulation with the observations in Figures 2, 6, 8 and 
9, it is evident that the redshift range z < 5.7 is consistent with being in the post-overlap stage 
of reionization, where Aj increases as the remaining overdense regions in the IGM are being 
ionized. The emergence of the GP trough at z ~ 6, and the stringent lower limits on the neutral 
hydrogen fraction from the presence of the Ly/3 and Lyy GP troughs suggest the onset of a rapid 
transition in the ionization state of the IGM. The upper limits on the ionizing background and 
mean free path, and the lower limits on the neutral fraction and critical overdensity we derive 
from the observations are all consistent with the typical values near the end of the overlap stage 
of reionization in the simulations. In this sense, the epoch of z ~ 6 is at the transition from the 
overlap stage to the post-overlap stage of reionization. Even though the current observations 
certainly cannot probe deeper into the reionization epoch, the near phase-transition behavior of 
the ionization state of the IGM at z ~ 6, and the narrow redshift range over which this process 
occurs in cosmological reionization simulations both imply that the reionization redshift cannot be 
at a redshift much higher than six. 

We caution, however, that these results are based on the GP trough in only one quasar, SDSS 
1030-1-0524. The IGM is unlikely to be reionized in an uniform manner. We have shown that at 
z ~ 6, the mean free path of the ionizing photons is likely to be shorter than 10 comoving Mpc. In 
this case, the ionizing background could show substantial fluctuations. These fluctuations could 
be due to the rarity and the possible strong clustering of ionizing sources, as well as the radiative 
transfer effects such as shadowing of the sources. For example, in the simulation of Gnedin (2000), 
the ionizing background shows fluctuation of a factor of a few at 2; ~ 6 (his Figure 5) in the 
underdensc regions of the IGM (which occupy most of the volume). In order to calculate the 
amplitude of these fluctuations, we would need detailed simulations of the formation of the first 
generation of ionizing objects. This is beyond the scope of the present paper The first generation 
of objects is likely to be highly biased and clustered, which would give rise to different reionization 
epochs along different lines of sight and in different regions of the IGM. If the universe was ionized 
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by luminous quasars (unlikely based on current statistical on quasar evolution, Fan et al. 2001c), 
the large HII regions produced by them near the epoch of overlap could even lead to gaps in 
transmitted flux (Haiman & Loeb 1998, Miralda-Escude et al. 2000). Thus, it is quite possible 
that the next quasar discovered at z > 6 will show somewhat different absorption properties. 



6. Discussion 

The presence of the GP trough in the spectrum of the highest redshift quasar provides a first 
opportunity to study structure formation at high-redshift and the formation of the first generation 
of stars and galaxies. To first order, the redshift of the reionization epoch depends on two factors: 
(a) the amount of small scale power in the power spectrum of mass density fluctuations, which 
determines how many halos can collapse to make stars at a certain redshift. It also determines 
the dumpiness of the IGM; and (b) how efficiently the collapsed halos can produce UV photons, 
which in turn depends on the stellar initial mass function, the efficiency of star formation, the 
escape factor of UV photons from proto-galaxies, etc. The reionization redshift and the observed 
ionizing background can be used to put constraints on these factors. 



In Chiu &: Ostriker (2000) 's semi-analytic calculations, the reionization redshift for low 
density A-dominated CDM models is typically in the range of 7 - 11, not very far from z ~ 6 — 8 
suggested in this paper. On the other hand, most of the non-Gaussian texture and isocurvature 
models have too much small scale power and a much earlier reionization redshift z ^ 10, making 
them difficult to reconcile with the observations presented in this paper. Barkana, Haiman & 
Ostriker (2001) used the constraint on the reionization epoch to put limits on models of Warm 
Dark Matter. Assuming that the reionization redshift is smaller than 7, we find that the mass of 
the warm dark matter particle must be smaller than 3 keV in the standard model they considered. 

Paper I estimates that the SDSS will find about 20 luminous quasars in the redshift range 
6 < z < 6.6 over the 10,000 deg^ survey area. High resolution, high signal-to-noise ratio spectra 
of these luminous quasars in the Lyman series absorption regions will provide valuable probes 
of the end of the reionization epoch, and in particular measure the spatial inhomogeneity of the 
reionization process. 

In order to probe deep into and beyond the reionization epoch, souces at higher redshift are 
needed. The Lya emission line moves out of the optical window and into the infrared at z ^ 6.6. 
Thus, the objects become very faint in the optical wavelengths due to Lyman series absorption by 
neutral hydrogen gas in the foreground at lower redshifts. Hence, this is the limiting redshift for 
discovering sources from ground-based optical surveys. Deep and wide-field IR surveys such as the 
PRIME0 mission are required to find objects at even higher redshifts. 



*http://prime. pha.jhu.edu/index. html 
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The stringent limits on the neutral hydrogen fraction at low redshifts from the GP test arise 
from the large cross section of neutral hydrogen to Lya photons. For the same reason, the GP test 
quickly becomes insensitive to larger neutral hydrogen fractions. With a half hour exposure on 
the Keck telescope, we are able to place a lower limit on the neutral fraction /hi ^ 10~^. As the 
neutral fraction hydrogen scales with r, or logarithmically with the flux, it is impossible to obtain 
a limit more than a few times stronger than the current limit for any reasonable exposure time. 
Therefore, the original version of the GP test cannot be used to probe deeper into the reionization 
epoch, for which /hi ~ 1. 

For an object observed prior to the reionization epoch, the damping wing of the GP trough 
arising from the large GP optical depth of the neutral medium will extend into the red side of 
the Lya emission line (Miralda-Escude 1998). The depth and extent of this damping wing can 
in principle be used to measure the neutral fraction and the redshift of the reionization epoch. 
However, this GP damping wing test cannot be applied to luminous quasars, due to the proximity 
effect from the quasar itself. As shown by Madau & Rees (2000) and Gen & Haiman (2000), the 
quasar will ionize a HIT region around it if the bulk of the IGM is still neutral, resulting in a line 
profile very similar to the case where the IGM has already been ionized. An alternative is to search 



for lower luminosity quasars (Stern et al. 2000) and galaxies ( Hu et al. 1999 , Rhoads Sz Malhotra 



2001 ) at high-redshifts by either applying the color dropout technique to deep multi-band imaging 
data or by search in regions around foreground galaxy clusters where background objects might 
be amplified by gravitational lensing (e.g Ellis et al. 2001| ). The highly magnified sources may be 



bright enough to allow high S/N spectroscopy, necessary for accurately measuring the profile of 
the damping wing. 

Another promising opportunity for probing the reionization epoch arises from the observations 
of IGM absorption in the afterglow of high redshift Gamma Ray Bursts (GRBs, Loeb 2001). 
GRBs are transient phenomena with no immediate impact on the surrounding IGM. The time 
dilation at high-redshift leads to a longer duration for the afterglows, making the identification and 
follow-up observations easier. Moreover, if high-redshift GRBs are associated with the collapse of 
massive stars, their evolution is likely to be similar to those of star-forming galaxies, which show 
slower decline at high redshift ( ^teidel et al. 1998| ) than the rapid decline in the number density 
of quasars (e.g. Fan et al. 2001b). 

We thank Renyue Gen, Zoltan Haiman, Avi Loeb, Pat McDonald, Jerry Ostriker, Martin 
Rees, Joop Schaye, and David Weinberg for helpful discussions. We acknowledge support from 
NSF grant PHYOO- 70928 and a Frank and Peggy TapUn Fellowship (XF), and NSF grant 
AST-0071091 (MAS). 
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Figure 1. Evolution of Lya absorption based on the observations of four quasars at 2; > 5.7 in 
Fan et al. (2001), Becker et al. (2001) and Pentericci et al. (2001). The results at ^abs < 5.6 
are averaged over four lines of sight, and the error bars include contributions from photon noise, 
uncertainty in the intrinsic quasar continuum, and an estimate of the sample variance. 
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Figure 2. Evolution of the photoionization rate (in units of 10^^^ s^^) with redshift. The points 
at 2;abs < 4.5 are taken from McDonald fc Miralda-Escude (2001) . For the highest redshift 
measurement (at z = 6.05), the upper limits based on Lya , Ly/3 and Ly7 GP troughs are shown 
separately. 
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Figure 3. Comparison of Lya forest spectra at (a) Zabs = 5.5, and (b) Zabs = 6.0. In each panel, the 
spectrum at the bottom shows the transmitted flux divided by a model for the continuum, in the 
Lya forest region at that redshift in the Keck spectrum of SDSS 1044-0125 and SDSS 1030+0524, 

respectively. The smaller chunks are simulated spectra along eight different lines of sight through 
the simulation cube, with the same resolution, binning and noise properties as the Keck spectra. 
Pairs of simulated spectra are offset by one unit in the vertical direction for clarity. 
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Figure 4. Probability distribution function (PDF) of the transmitted flux in the Lya forest region 
in the Keck spectra (points) and simulated spectra (lines) at (a) Zabs = 5.5, (b) Zabs = 5.7 and (c) 
Zabs = 6.0. The solid line shows the average flux PDF computed using 400 artificial spectra, and 
the dashed lines shows the expected 1 — a uncertainty in the PDF from a single spectrum of length 
equal to that of the real spectrum (i.e., the cosmic variance). 
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Figure 5. Threshold crossing statistic as a function of the fraction of pixels in the spectrum less 
than the flux threshold in the Lya forest region in the Keck spectra (points) and simulated spectra 
(lines) at (a) Zabs = 5.5, (b) Zabs = 5.7 and (c) z^hs = 6.0. The solid line shows this statistic 
computed using 400 artificial spectra, and the dashed lines shows the expected 1 — a uncertainty in 
this quantity from a single spectrum of length equal to that of the real spectrum (i.e., the cosmic 
variance) . 
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Figure 6. Evolution of neutral hydrogen fraction of the IGM. The solid points with error 
bars are measurements based on the four high-redshift quasars, for both mass-averaged, and 

volume-averaged neutral fraction. The solid and dashed lines are the mass- and volumc-avcragcd 
results from the N128_L4_A simulation of Gnedin (2000). The neutral fraction inferred from the 
observations is comparable to that of the transition from overlapping stage to post-overlap stage 
of reionization in the simulation. 
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Figure 7. Cumulative distribution of transmitted flux ratio as a function of density at z = 5.3, 
with a transmitted flux ratio of 0.08, corresponding to r_i2 = 0.15. The dashed line shows the 
probability distribution function of the density. Most of the transmitted flux comes from the 
under dense voids in the IGM. 
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Figure 8. Evolution of the limiting density, Aj, defined as the density above which 95% of the HI 
gas lies. In this picture, reionization starts in voids, and gradually penetrates deeper into overdense 
regions in an inhomogeneous IGM. Both results inferred from the observations (solid line) and those 
from the simulation of Gnedin (2000, dashed line) are shown. 
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Figure 9. Evolution of the mean free path of ionizing photons inferred from the observations (sohd 
hne) and simulation of Gnedin (2000, dashed line). 



